#ifndef MATRIX_H
#define MATRIX_H

#include "rect.h"

#include <assert.h>
#include <iostream>
#include <vector>
#include <malloc.h>
#include <memory.h>

/*!
 * \brief The Matrix class contain a set of members and methods to manipulate a two dimentional matrix
 */
template <typename T>
class Matrix
{
    T* _data;

public:

    int _height;/*!< \brief represent the matrix row count */
    int _width; /*!< \brief represent the matrix column count */

    explicit Matrix() :
        _height(0),
        _width(0),
        _data(NULL)
    {
    }

    /*! \brief Matrix copie constructor*/
    Matrix(const Matrix& value) :
        _height(0),
        _width(0),
        _data(NULL)
    {
        (*this) = value;
    }

    /*!
     * \brief Matrix instanciate a matrix of with `height` colmuns and `width` rows.
     * \param width matrix column count
     * \param height matrix row count
     */
    explicit Matrix(unsigned int width,
                    unsigned int height) :
        _width(width),
        _height(height)
    {
        _data = new T[_width * _height]();
    }


    /*!
     * \brief Matrix instanciate a matrix of with `height` colmuns and `width` rows.
     * \param width matrix column count
     * \param height matrix row count
     * \param initial value to initializate matrix indices with
     */
    explicit Matrix(unsigned int width,
                    unsigned int height,
                    T initial) :
        _width(width),
        _height(height)
    {
        _data = new T[_width * _height]();
        for(unsigned int i=0; i< _width; i++)
            for(unsigned int j=0; j< _height; j++)
                *at(i,j) = initial;
    }

    /*!
     * \brief Matrix instanciate a matrix of `height` colmuns and `width` rows.
     * \param width matrix column count
     * \param height matrix row count
     * \param data pointer to the matrix contained data (of type `T`). Sould be an instanciated array of size `width * height`
     */
    explicit Matrix(unsigned int width,
                    unsigned int height,
                    const T* data) :
        _width(width),
        _height(height)
    {
        assert(_data != ((void*)0));
        _data = data;
    }

    /*!
     * \brief Matrix instanciate a square matrix with a size of `diagonal`.
     * \param diagonal strait line matrix size
     * \param initial value to initializate matrix indices with
     */
    explicit Matrix(unsigned int diagonal,
                    T initial) :
        _width(diagonal),
        _height(diagonal)
    {
        assert(diagonal> 0);
        _data = new T[diagonal * diagonal]();
        for(unsigned int i=0; i< _width; i++)
            for(unsigned int j=0; j< _height; j++)
                *at(i,j) = initial;
    }

    /*!
     * \brief Matrix instanciate a square matrix with a size of `diagonal`.
     * \param diagonal strait line matrix size
     * \param data pointer to the matrix contained data (of type `T`). Sould be an instanciated array of size `diagonal * diagonal`
     */
    explicit Matrix(unsigned int diagonal,
                    const T* data) :
        _width(diagonal),
        _height(diagonal)
    {
        assert(diagonal> 0);
        assert(_data != ((void*)0));
        _data = data;
    }

    ~Matrix()
    {
        delete _data;
    }

    /*! \brief at return a mutable value of the matrix at `row` and `col`*/
    inline T* at (unsigned int row,
           unsigned int col) {
        return &_data[row * _width + col];
    }

    /*! \brief at return a constant value of the matrix at `row` and `col`*/
    inline const T* at(unsigned int row,
                unsigned int col) const {
        return &_data[row * _width + col];
    }

    inline static void multiplyAdd(const Matrix& a, const Matrix& b, Matrix& c) {
        assert(a.width() == b.height());
        assert((c.width() == b.width()) && (c.height() == a.height()));
        //NOTE: could be optimized with SIMD instruction with template class specialization. (BLAS, Lapack)
        //NOTE: could be parallelized using openMP.
        for(unsigned int i=0; i < b.width(); i++)
            for(unsigned int j=0; j < a.height(); j++)
                for(unsigned int k=0; k<a.width(); k++)
                    *c.at(i,j)+= *a.at(i,k) * *b.at(k,j);
    }

    /*! \brief operator = Matrix affection operator */
    inline Matrix& operator=(const Matrix& value)
    {
        if(&value != this) {
            if(_data == NULL) {
                _data = new T[value.width() * value.height()]();
                _height = value.height();
                _width = value.width();
            }
            else if((_width != value.width()) || (_height != value.height())) {
                _data = (T*)realloc((void*)_data,sizeof(T) * value.height() * value.width());
                _height = value.height();
                _width = value.width();
            }
            memcpy((void*)_data,(void*)value.data(),sizeof(T) * value.width() * value.height());
        }
        return *this;
    }

    /*!  \brief operator () return the current matrix part contained in the rectangle `rect` */
    inline Matrix operator()(const Rect& rect) const
    {
        assert((rect.w + rect.x) < _width);
        assert((rect.h + rect.y) < _height);
        assert((rect.h > 0) && (rect.w > 0));
        Matrix ret(rect.w,rect.h);
        //NOTE: could be optimized with SIMD instruction with template class specialization.
        for (int i = 0; i < rect.w; ++i)
            for (int j = 0; j < rect.h; ++j)
                *ret.at(i,j) = *at(i + rect.x, j + rect.y);
        return ret;
    }
    /*!  \brief operator () return the current matrix part contained in the rectangle from `(0,0)` to `(width,height)` */
    inline Matrix operator()(unsigned int width,
                             unsigned int height) const{
        Rect rect;
        rect.x = 0;
        rect.y = 0;
        rect.w = width;
        rect.h = height;
        return (*this)(rect);
    }

    /*!
     * \brief operator * return the matrix product between current matrix and `value` matrix.
     * The `value` matrix height sould be equal to the current matrix width.
     */
    inline Matrix operator* (const Matrix& value) const
    {
        assert(width() == value.height());
        Matrix mat(value.width(),height());
        //NOTE: could be optimized with SIMD instruction with template class specialization. (BLAS, Lapack)
        //NOTE: could be parallelized using openMP.
        for(unsigned int i=0; i<value.width(); i++) {
            for(unsigned int j=0; j<height(); j++) {
                *mat.at(i,j) = *at(i,0) * *value.at(0,j);
                for(unsigned int k=1; k<_width; k++) {
                    *mat.at(i,j)+= *at(i,k) * *value.at(k,j);
                }
            }
        }
        return mat;
    }

    /*!
     * \brief operator * return the matrix product between current matrix and `value` matrix.
     * The `value` matrix height sould be equal to the current matrix width.
     */
    inline Matrix& operator*= (const Matrix& value)
    {
        *this = *this * value;
        return *this;
    }

    /*!
     * \brief operator + add each indices of current matrix with the `value` matrix.
     * The `value` should have same dimension than the current matrix.
     */
    inline Matrix operator+(const Matrix& value) const
    {
        assert((width() == value.width()) && (height() == value.height()));
        Matrix mat(width(),height());
        for(unsigned int i=0; i<value.width(); i++)
            for(unsigned int j=0; j<height(); j++)
                *mat.at(i,j) = *at(i,j) + *value.at(i,j);
        return mat;
    }

    /*!
     * \brief operator - substract each indices from the current matrix with the `value` matrix.
     * The `value` should have same dimension than the current matrix.
     */
    inline Matrix operator-(const Matrix& value) const
    {
        assert((width() == value.width()) && (height() == value.height()));
        Matrix mat(width(),height());
        for(unsigned int i=0; i<value.width(); i++)
            for(unsigned int j=0; j<height(); j++)
                *mat.at(i,j) = *at(i,j) - *value.at(i,j);
        return mat;
    }

    /*! \brief operator << write the matrix in the stream */
    inline friend std::ostream& operator<<( std::ostream &stream,
                                            const Matrix& value)
    {
        stream << value.width() << " " << value.height() << std::endl;
        for(unsigned int i=0; i<value.width(); i++) {
            for(unsigned int j=0; j < (value.height() - 1); j++)
                stream << *value.at(i,j) << ";";
            stream << *value.at(i,value.height() - 1) << std::endl;
        }
        return stream;
    }

    /*! \brief operator >> instanciate the matrix trought the stream */
    inline friend std::istream& operator>>( std::istream &stream,
                                            Matrix& value)
    {
        unsigned int w, h;
        stream >> w >> h;
        value = Matrix<T>(w,h);
        T* ptr = value._data;
        while(stream.peek() !=  std::char_traits<char>::eof()) {
            stream >> *ptr;
            ptr++;
            stream.get();
        }
        return stream;
    }

    /*! \brief height return the matrix height */
    inline int height() const {
        return _height;
    }

    /*! \brief width return the matrix width */
    inline int width() const {
        return _width;
    }

    /*! \brief return the data matrix pointer */
    inline  T* data()  {
        return _data;
    }

    /*! \brief return the data matrix pointer */
    inline  const T* data() const {
        return _data;
    }


};

#endif // MATRIX_H
